home *** CD-ROM | disk | FTP | other *** search
/ The Arsenal Files 4 / The Arsenal Files 4 (Arsenal Computer).ISO / ham / sattrk31.tgz / sattrack-3.1.tar / SatTrack / src / sattrack / satrefr.c < prev    next >
C/C++ Source or Header  |  1995-03-16  |  17KB  |  298 lines
  1. /******************************************************************************/
  2. /*                                                                            */
  3. /*  Title       : satrefr.c                                                   */
  4. /*  Author      : Manfred Bester                                              */
  5. /*  Date        : 01Jan92                                                     */
  6. /*  Last change : 01Jan95                                                     */
  7. /*                                                                            */
  8. /*  Synopsis    : This function calculates the correction to elevation due    */
  9. /*                to atmospheric refraction for the wavelength region between */
  10. /*                200 and 20000 nm.                                           */
  11. /*                                                                            */
  12. /*                                                                            */
  13. /*  Input parameters :                                                        */
  14. /*                                                                            */
  15. /*   1.  standard atmospheric pressure      [mmHg]   :  atmPress              */
  16. /*   2.  standard ambient temperature       [deg C]  :  ambTemp               */
  17. /*   3.  relative humidity                  [%]      :  relHumid              */
  18. /*   4.  wavelength (in vacuum)             [micron] :  waveLength            */
  19. /*   5.  reference wavelength (in vacuum)   [micron] :  refWaveLength         */
  20. /*   5.  Earth's effective radius           [km]     :  effEarthRad           */
  21. /*   6.  scale height of the atmosphere     [km]     :  scaleHeight           */
  22. /*   7.  altitude of the observer           [km]     :  obsAlt                */
  23. /*   8.  elevation of the satellite         [rad]    :  satEle                */
  24. /*                                                                            */
  25. /*                                                                            */
  26. /*  Output parameters :                                                       */
  27. /*                                                                            */
  28. /*   1.  index of refraction                [---]    :  refIndex              */
  29. /*   2.  correction to elevation            [rad]    :  refrEle               */
  30. /*   3.  differential elevation correction  [rad]    :  diffEle               */
  31. /*   4.  air mass                           [---]    :  airMass               */
  32. /*                                                                            */
  33. /*                                                                            */
  34. /*  For reference see :   C.W.Allen, "Astrophysical Quantities",              */
  35. /*                        The Athlone Press, London (1976), pages 122, 124    */
  36. /*                                                                            */
  37. /*                        P.N.Tverskoi, "Physics of the Atmosphere, A Course  */
  38. /*                        in Meteorology", E.S.Selezneva (ed.), GIMIZ,        */
  39. /*                        Leningrad (1962)                                    */
  40. /*                        (for conversion of relative humidity into water     */
  41. /*                        vapor pressure)                                     */
  42. /*                                                                            */
  43. /*                                                                            */
  44. /*  Original papers : a)  H.Barrel, J. Opt. Soc. Am. 41, 295-299 (1951)       */
  45. /*                    b)  B.Edlen,  J. Opt. Soc. Am. 43, 339-344 (1953)       */
  46. /*                    c)  C.D.Coleman, W.R.Bozman, W.F.Meggers, "Table of     */
  47. /*                        Wavenumbers", NBS Monograph 3, Washington D.C.      */
  48. /*                        (1960)                                              */
  49. /*                                                                            */
  50. /*  All equations were taken from the original papers :                       */
  51. /*                                                                            */
  52. /*  assumed standard temperature and pressure: +15 deg C, 760 mmHg            */
  53. /*                                                                            */
  54. /*  n1    = (n(15,760)-1)*1.0e-6          from eqn. (2) in b)                 */
  55. /*  wv    = correction for humid air      from eqn. (5) in b)                 */
  56. /*  ntpf1 = (n(t,p,f) -1)*1.0e-6          from eqn. (9) in a)                 */
  57. /*                                                                            */
  58. /*  reference (c) has been used for cross-checking the formulae and the data  */
  59. /*                                                                            */
  60. /*                                                                            */
  61. /*  SatTrack is Copyright (c) 1992, 1993, 1994, 1995 by Manfred Bester.       */
  62. /*  All Rights Reserved.                                                      */
  63. /*                                                                            */
  64. /*  Permission to use, copy, and distribute SatTrack and its documentation    */
  65. /*  in its entirety for educational, research and non-profit purposes,        */
  66. /*  without fee, and without a written agreement is hereby granted, provided  */
  67. /*  that the above copyright notice and the following three paragraphs appear */
  68. /*  in all copies. SatTrack may be modified for personal purposes, but        */
  69. /*  modified versions may NOT be distributed without prior consent of the     */
  70. /*  author.                                                                   */
  71. /*                                                                            */
  72. /*  Permission to incorporate this software into commercial products may be   */
  73. /*  obtained from the author, Dr. Manfred Bester, 1636 M. L. King Jr. Way,    */
  74. /*  Berkeley, CA 94709, USA. Note that distributing SatTrack 'bundled' in     */
  75. /*  with ANY product is considered to be a 'commercial purpose'.              */
  76. /*                                                                            */
  77. /*  IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, */
  78. /*  SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF   */
  79. /*  THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN ADVISED  */
  80. /*  OF THE POSSIBILITY OF SUCH DAMAGE.                                        */
  81. /*                                                                            */
  82. /*  THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT      */
  83. /*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A   */
  84. /*  PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"      */
  85. /*  BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, */
  86. /*  UPDATES, ENHANCEMENTS, OR MODIFICATIONS.                                  */
  87. /*                                                                            */
  88. /******************************************************************************/
  89.  
  90. #include <stdio.h>
  91. #include <math.h>
  92.  
  93. #ifndef STDLIB
  94. #include <stdlib.h>
  95. #endif
  96.  
  97. #include "satglobalsx.h"
  98. #include "sattrack.h"
  99.  
  100. /******************************************************************************/
  101. /*                                                                            */
  102. /* getRefraction: calculates the refraction in the atmosphere                 */
  103. /*                                                                            */
  104. /******************************************************************************/
  105.  
  106. void getRefraction(atmPress,ambTemp,relHumid,waveLength,refWaveLength,
  107.                    effEarthRad,scaleHeight,obsAlt,satEle,
  108.                    prefIndex,prefrEle,pdiffEle,pairMass)
  109.  
  110. double atmPress, ambTemp, relHumid, waveLength, refWaveLength, effEarthRad; 
  111. double scaleHeight, obsAlt, satEle, *prefIndex, *prefrEle, *pdiffEle, *pairMass;
  112.  
  113. {
  114.     double ntpf1, ntpf1r, n1, n1r, ntpf, ntpfr, mr, mw, nr, tp, wv, wvr;
  115.     double waterVapor, absTemp, log10E, E;
  116.     double refrEle, diffEle, airMass, hAtm, sinEle, arg;
  117.     double rSin, rSin2, rObs, rObs2, rAtm, rAtm2, a, alpha, sinAlp, beta, betar;
  118.  
  119. /******************************************************************************/
  120. /*                                                                            */
  121. /* return zero if satellite is below the horizon                              */
  122. /*                                                                            */
  123. /******************************************************************************/
  124.  
  125.     if (satEle < 0.0)
  126.     {
  127.         *prefIndex = 0.0;
  128.         *prefrEle  = 0.0;
  129.         *pdiffEle  = 0.0;
  130.         *pairMass  = 0.0;
  131.  
  132.         return;
  133.     }
  134.  
  135. /******************************************************************************/
  136. /*                                                                            */
  137. /* convert relative humidity [%] into water vapor pressure [mmHg]             */
  138. /*                                                                            */
  139. /* scale factor : 1/133.322 = 1/100/1.33322                                   */
  140. /* (conversion from "per cent" into decimal number and from mbar into mmHg)   */
  141. /*                                                                            */
  142. /* for reference see :   P.N.Tverskoi, "Physics of the Atmosphere, A Course   */
  143. /*                       in Meteorology", E.S.Selezneva (ed.), GIMIZ,         */
  144. /*                       Leningrad (1962)                                     */
  145. /*                                                                            */
  146. /******************************************************************************/
  147.  
  148.     absTemp    = 273.155 + ambTemp;
  149.     log10E     = 23.5518 - 2937.4 / absTemp - 4.9283 * log10(absTemp); 
  150.     E          = pow(10.0,log10E);
  151.     waterVapor = E * relHumid / 133.322;
  152.  
  153. /******************************************************************************/
  154. /*                                                                            */
  155. /* rescale length units from km into m                                        */
  156. /*                                                                            */
  157. /******************************************************************************/
  158.  
  159.     effEarthRad *= CKMM;
  160.     scaleHeight *= CKMM;
  161.     obsAlt      *= CKMM;
  162.  
  163. /******************************************************************************/
  164. /*                                                                            */
  165. /* first calculate (n-1) for standard temperature (t = +15 deg C) and         */
  166. /* standard pressure (p = 760 mmHg) in dry air for both the wavelength and    */
  167. /* the reference wavelength                                                   */
  168. /*                                                                            */
  169. /* waveLength and refWaveLength are the two wavelengths in vacuum             */
  170. /*                                                                            */
  171. /******************************************************************************/
  172.  
  173.     mw  = 1.0 / (waveLength * waveLength);
  174.     mr  = 1.0 / (refWaveLength * refWaveLength);
  175.  
  176.     n1  = 64.328 + 29498.1 / (146.0 - mw) + 255.4 / (41.0 - mw);
  177.     n1r = 64.328 + 29498.1 / (146.0 - mr) + 255.4 / (41.0 - mr);
  178.  
  179. /******************************************************************************/
  180. /*                                                                            */
  181. /* now calculate the contribution due to water vapor in the air for both      */
  182. /* wavelengths                                                                */
  183. /*                                                                            */
  184. /******************************************************************************/
  185.  
  186.     tp  = 1.0 + 0.003661 * ambTemp;
  187.  
  188.     wv  = -(0.0624 - 0.000680 * mw) * waterVapor / tp;
  189.     wvr = -(0.0624 - 0.000680 * mr) * waterVapor / tp;
  190.  
  191. /******************************************************************************/
  192. /*                                                                            */
  193. /* calculate pressure at the specified altitude                               */
  194. /*                                                                            */
  195. /* then calculate (n-1) for pressure (atmPress) and temperature (ambTemp)     */
  196. /*                                                                            */
  197. /* explanation of the numerical factor :                                      */
  198. /*                                                                            */
  199. /* 1.387188936e-3 = (1 + 15*0.003661)/(760*(1 + 760*0.813e-6))                */
  200. /*                                                                            */
  201. /******************************************************************************/
  202.  
  203.     nr     = (1.0 + (1.049 - 0.0157 * ambTemp) * 1.0e-6 * atmPress) / tp;
  204.     nr    *= 1.387188936e-3;
  205.  
  206.     ntpf1  = n1  * atmPress * nr + wv;
  207.     ntpf1r = n1r * atmPress * nr + wvr;
  208.  
  209.     ntpf   = 1.0 + ntpf1  * 1.0e-6;
  210.     ntpfr  = 1.0 + ntpf1r * 1.0e-6;
  211.  
  212. /******************************************************************************/
  213. /*                                                                            */
  214. /* calculate the angle of normal incidence into the atmosphere (alpha)        */
  215. /* and then the correction to elevation (refralt) by applying Snell's law:    */
  216. /*                                                                            */
  217. /* n2/n1 = ntpf = sin(alpha)/sin(beta)  and  refrEle = alpha - beta           */
  218. /*                                                                            */
  219. /* if the elevation is still or becomes less than 0.0 deg the correction      */
  220. /* (refrEle) is set to zero                                                   */
  221. /*                                                                            */
  222. /******************************************************************************/
  223.  
  224.     rObs    = effEarthRad + obsAlt;
  225.     rAtm    = effEarthRad + scaleHeight;
  226.     hAtm    = scaleHeight - obsAlt;
  227.  
  228.     sinEle  = sin(satEle);
  229.     rSin    = rObs * sinEle;
  230.     rSin2   = rSin * rSin;
  231.     rObs2   = rObs * rObs;
  232.     rAtm2   = rAtm * rAtm;
  233.  
  234.     a       = -rSin + sqrt (rSin2 + rAtm2 - rObs2);
  235.     airMass = a / hAtm;
  236.  
  237.     alpha   = acos((a*a + rAtm2 - rObs2) / (2.0*a*rAtm));
  238.     sinAlp  = sin(alpha);
  239.  
  240.     arg     = sinAlp / ntpf;
  241.  
  242.     if (arg >  1.0) arg =  1.0;      /* protect against out-of-range problems */
  243.     if (arg < -1.0) arg = -1.0;
  244.  
  245.     beta    = asin(arg);
  246.     arg     = sinAlp / ntpfr;
  247.  
  248.     if (arg >  1.0) arg =  1.0;
  249.     if (arg < -1.0) arg = -1.0;
  250.  
  251.     betar   = asin(arg);
  252.  
  253.     refrEle = alpha - beta;
  254.     diffEle = beta  - betar;
  255.  
  256. /******************************************************************************/
  257. /*                                                                            */
  258. /* correction for refraction in the atmosphere is now done with respect to:   */
  259. /*                                                                            */
  260. /* 1. wavelength of observation                                               */
  261. /* 2. ambient temperature                                                     */
  262. /* 3. atmospheric pressure                                                    */
  263. /* 4. water vapor pressure                                                    */
  264. /* 5. altitude (zenith distance)                                              */
  265. /* 6. curvature of the atmosphere                                             */
  266. /*                                                                            */
  267. /******************************************************************************/
  268.  
  269. /*
  270.     printf("\n");
  271.     printf("wavelength              [micron] : %9.6f\n",waveLength);
  272.     printf("reference wavelength    [micron] : %9.6f\n",refWaveLength);
  273.     printf("atmospheric pressure    [mmHg]   : %9.3f\n",atmPress);
  274.     printf("ambient temperature     [deg C]  : %9.3f\n",ambTemp);
  275.     printf("water vapor pressure    [mmHg]   : %9.3f\n",waterVapor);
  276.     printf("refractive index        (n-1)e+6 : %9.3f\n",ntpf1);
  277.     printf("alpha                   [deg]    : %9.6f\n",alpha*CRD);
  278.     printf("beta                    [deg]    : %9.6f\n",beta*CRD);
  279.     printf("betar                   [deg]    : %9.6f\n",betar*CRD);
  280.     printf("refractive correction   [arcsec] : %9.4f\n",refrEle*CRAS);
  281.     printf("differential correction [arcsec] : %9.4f\n",diffEle*CRAS);
  282.     printf("air mass                [---]    : %9.3f\n",airMass);
  283. */
  284.  
  285.     *prefIndex = ntpf; 
  286.     *prefrEle  = refrEle;
  287.     *pdiffEle  = diffEle;
  288.     *pairMass  = airMass;
  289.  
  290.     return;
  291. }
  292.  
  293. /******************************************************************************/
  294. /*                                                                            */
  295. /* End of function block satrefr.c                                            */
  296. /*                                                                            */
  297. /******************************************************************************/
  298.